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ABSTBACT 

A  new  method  is  described  for  computing  the  effect  of  correlation,  in- 
homogeneity  and  exchange  on  the  Thomas-Fermi  model  of  the  atom.   The  method 
makes  use  of  the  many  body  point  of  view,  rather  than  an  independent  parti- 
cle point  of  vJ ew,  by  considering  the  hierarchy  equation  linking  the  n  par- 
ticle Green's  functions.   The  hierarchy  is  truncated  by  a  prescription  equi- 
valent to  the  Gell-Mann  and  Brueckner  theory  of  the  high  density  electron 
gas,  resulting  in  a  description  of  the  atom  in  which  the  exchange  interac- 
tion is  replaced  by  the  effective  interaction.   The  physical  significance 
of  this  replacement  is  noted. 

The  Green's  function  for  this  model  is  then  expanded  as  a  series  in 
powers  of  ^.      The  lowest  order  term  is  found  to  describe  the  Thomas -Fermi 
model  of  the  atom.   The  equation  for  the  next  higher  term  contributing  to 
this  expansion  is  manipulated  so  as  to  yield  an  ordinary  differential  equa- 
tion for  the  corresponding  correction  to  the  potential.   This  equation  con- 
tains a  term  which  expresses  the- effect  of  inhomogeneity  and  another  which 
arises  from  the  correlation  of  the  electrons  and  from  exchange.   The  inhomo- 
geneity term  is  one  which  has  been  found  previously.   Study  of  the  correla- 
tion term  shows  that  it  depends  on  the  separation  energy  of  an  electron  from 
an  infinite  electron  gas,  which  suggests  a  generalization  by  which  the  me- 
thod might  be  made  applicable  to  those  outer  regions  of  the  atom  for  which 
the  electron  density  is  below  that  to  which  the  Gell-Mann  and  Brueckner 
theory  would  apply. 
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INTRODUCTION 


In  a  previous  paper,   we  presented  a  systematic  method  for  deriving  the 
Thom.as-Fermi  equation  and  quantum  corrections  (exchange  andinhomogeneity) 
from  the  many  body  equations  of  motion  of  an  atomic  system.   This  method  vas 
based  on  the  Green's  function  formulation  of  the  many  body  problem,   and,  in 
particular,  made  use  of  the  Hartree-Fock  approximation  vlthin  that  formalism. 
We  shoved  that  the  single  particle  Green's  function  in  the  Hartree-Fock 
approximation  could  be  expanded  as  a  series  in  powers  of  ^,  and  that  the  low- 
est order  term  in  this  series  described  the  Thomas-Fermi  distribution.   The 
next  non-vanishing  higher  order  term  was  found  to  be  the  same  as  the  quantum 
correction  foimd  by  Kompaneets  and  Pavlovskii,   who  had  based  their  work  on 
the  Hartree-Fock  approximation  in  the  density  matrix  formalism,  and  indcpen- 
dently  by  Kirzhnits,   who  had  used  a  method  involving  successive  commutators, 
also  within  the  Hartree-Fock  scheme. 

The  effect,  however,  of  including  all  the  higher  order  terms  of  our  in- 
finite series  can  be  to  recover  only  the  Hartree-Fock  solution,  and  conse- 
quently neglects  dynamic  particle  correlation.   In  this  paper,  we  shall  use 
the  technique  developed  in  I  to  include  some  of  the  correlation  effects  which 
are  discarded  ab  Initio  in  a  Hartree-Fock  approach.   We  shall  do  this  by  us- 
ing  the  Gell-Mann-Brueckner  theory  of  the  high  density  electron  gas^  as  the 
basic  approximation  instead  of  the  Hartree-Fock  approximation.   It  will  still 
be  possible  to  express  the  single  particle  Green's  function  as  a  series  in 
powers  of -fi,  and  Just  as  before,  the  lowest  order  term  will  describe  the 
Thomas-Fermi  distribution.   This  time,  however,  the  higher  order  terms  which 
may  be  regarded  as  corrections  to  the  Thomas-Fermi  model  introduce  not  only 


the  quantum  corrections  but  also  those  correlation  effects  vith  which  Gell- 
Mann  and  Brueckner  were  concerned. 

The  logical  approach  to  this  problem  is  divided  into  four  parts^  each 
of  which  is  the  subject  of  one  of  the  following  sections  of  this  paper.   In 
Section  II,  we  formulate  the  many  body  problem  in  the  mixed  position-momen- 
tum representation  of  the  Green's  function.   The  description  of  the  interac- 
tion between  particles  is  contained  formally  in  the  irreducible  self  energy 
operator.    We  choose  an  approximate  form  for  this  operator  which  describes 
the  important  interactions  of  the  G.B.  theory  and  find  that  the  only  differ- 
ence between  this  approximation  and  that  which  corresponds  to  the  Hartree- 
Fock  method  is  the  appearance  of  a  time  dependent  effective  interaction 
instead  of  the  instantaneous  coulomb  potential  in  the  exchange  term.   The 
direct  interaction  term  is  unaltered.   This  result  has  physical  signifi- 
cance which  will  guide  our  later  steps. 

In  Section  III,  we  expand  the  quantities  of  interest  as  series  in 
powers  of  if,  employing  the  techniques  presented  in  I.   By  series  in  powers 
of  xh",  we  mean  that  as  ^  is  considered  to  be  made  smaller,  the  lowest  order 
terms  dominate  higher  order  terms  to  greater  and  greater  extent.   Somce  care 
is  required  in  expanding  the  effective  interaction  in  order  to  avoid  recov- 
ering a  series  similar  to  the  one  summed  by  G.B.   We  can  obtain  a  useful 
expansion  of  the  effective  interaction  if  we  recognize  why  the  straightfor- 
ward expansion  fails  to  work  and  proceed  accordingly.   The  expansion  which 
results  is  completely  in  accord  with  the  physical  significance  of  the  result 
we  noted  in  the  preceding  section. 

In  Section  IV,  we  obtain  equations  for  the  lowest  order  terms  in  the 
expansion.  We  find  that  the  zeroth  order  term  describes  the  Thomas-Fermi 
distribution  and  that  the  first  order  correction  term  vanishes,  just  as  in  I. 


The  equation  for  the  second  order  term  is  quite  similar  to  that  found  in  I, 
the  sole  difference  being  in  that  term  which  in  I  arose  from  the  exchange 
term  and  here  arises  from  the  effective  interaction.   We  study  this  term  in 
detail  in  Section  V  and  find  that  it  represents  the  non-kinetic  part  of  the 
separation  energy  in  both  cases.   This  suggests  a  way  to  correct  our  work 
for  an  approximation  made  earlier  and,  also,  for  the  approximation  made  in 
using  the  G. B.  correlation  corrections,  which  are  valid  only  in  the  high  den- 
sity limit. 

Section  VI  contains  a  siMimary  and  a  discussion  of  the  similarities  and 
differences  between  this  work  and  an  earlier  paper  by  Lewis   on  essentially 
the  same  problem.   Our  result  in  this  work,  then,  is  an  equation  which  deter- 
mines the  lowest  correction  to  the  Thomas-Fermi  model  for  the  effect  of  cor- 
relation, exchange  and  inhomogeneity.   Solutions  to  this  equation  for  various 
atoms  will  be  deferred  to  a  future  publication. 

II  FORMULATION  OF  THE  PROBLEM 
We  again  consider  an  atom  or  an  ion  to  consist  of  a  nucleus  with  charge 
Ze  and  infinite  mass  located  at  the  origin,  surrounded  by  N  electrons,  each  of 
mass  m.   We  neglect  spin  dependent  forces  and  all  relativistic  effects  in  our 
system.   For  such  a  system,  the  equation  relating  the  one  particle  ajid  two 
particle  Green's  function  is 


i^(a/St^)  -  H^(r^)jG^(l;l') 

+  i  rd(2)v(l-2)G2(l,2;l'2"^)  =/65(l-l')  (2.1a) 


H  (r)  =  (-fi^/2m)V^  -  Ze^/lr|  (2.1b) 
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v(l-2)  =  5(t^-t2)e^/|r^-r2|.  (2.1c) 


This  notation  is  the  same  as  in  I:  each  numerical  argument  1,  1',  etc.  speci- 
fies a  spatial  coordinate,  a  time, and  a  spin  index.   Integration  over  a  numer- 
ical argument  implies  summation  over  the  associated  spin  index;  the  delta 
function  of  numerical  indices  signifies  the  product  of  a  Dirac  delta  function 
for  each  of  the  space  and  time  coordinates  with  a  Kronecker  delta  function  of 
the  spin  index.   The  symbol  +  as  a  superscript  to  a  numerical  index  indicates 
that  the  time  associated  vith  that  index  is  to  be  made  later  by  an  infinitesi- 
mal amount. 

The  one  and  tvo  particle  Green's  functions  G^  and  Gp  are  those  defined  in 
I  (c.f.  I,  5-7) •   They  are  of  use  in  computing  the  ground  state  expectation 
value  of  any  one  particle  or  tvo  particle  operator.   For  example,  the  density 

of  particles  at  spatial  point  r  having  spin  index  a,  is  given  by 

.v+i  1 

n(l)  =  -iG^Cljl"*"). 

The  ground  state  energy  E  is  also  easily  determined  (c.f.  I,  2.5). 

6 
Klein  and  Prange  have  introduced  the  irreducible  self  energy  operator  M 

by  giving  an  equation  defining  its  effect: 

i  /'d(2)v(l-2)G2(l.2;l',2+)  =  -  f  d(2)M(l,2)G(2;l').  (2.2) 

We  shall  use  this  definition  to  define  a  kernel  K(l,2), 

K(l,2)  =   i*(a/at^)  -  H^(r^)  6(l-2)-M(l,2)  ''       ■■  (2.5) 

such  that  Eq.  (2.1)  takes  the  form 

rd(2)K(l,2)G^(2;l')  =*5(1-1').   "' ■■   "'     ^    '"'".'■     '  (2.4) 
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spin 
For  a  system  such  as  ours  in  which^dependent  forces  are  neglected;,  a 

particle  propagates  without  change  of  spin.   That  is,  the  single  particle 

Green's  function  vanishes  unless  its  two  spin  indices  are  equal,  but  it  is 

otherwise  independent  of  the  spin.   Hence,  we  need  consider  Eq.  (2.i|)  only 

for  the  case  a,  =  cTp  =  a,  ,  .   Let  us  transform  this  equation  to  the  mixed 

position-momentum  representation  by  taking  the  Fourier  transform  with  respect 

to  r  -r  '  and  t  -t  '  ,holding  R  =  ^(r  +r  ' )  fixed,  denoting  the  transform 

variable's  by  p  and  w. 


J  d(r^-r^')d(t^-t^')d(2)K(l,2)G^(2,l' 


x^^H-HS'^rrii')  -"^V\' 


/-h>  =  h. 


This  transformation  was  used  extensively  in  I.  where  it  was  shown  that  the 
left  side  of  the  preceding  equation  could  be  expressed  in  term^s  of  the  trans- 
forms of  K  and  G  by  means  of  a  differential  operator  of  infinite  order,  0  as 


e 


K(R,pu),G(R,pu)) 


The  transforms  K  and  G  are  defined  by 


K(R,pu))    =   /    d(r  -r    )d(t   -t    )K(1,2)      _ 
«,  w  J        "~x  ■^<i         X     d.  a, -Op 


X  exp  < -1 


p.(r^-r^)    -   a)(t   -t    ) 

A»A  **'  art   ^  X  <— 


/fi 


(2.5) 


G(R,pa,)    =        d(„Vr2)d(t^-t2)G^(l,2)      ^^ 


1~"2 


X  exp  ^-l|p.(r  -r^)    -  u)(t   -t    )    /ft 


(2.6) 
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(The  operator  ©  is  defined  fully  and  its  effect  given  explicitly  in  I.   Its 
first  fev  terms  are  given  in  Appendix  B  of  this  paper.)   In  terms  of  this 
notation,  the  transform  of  Eq.  {2.k)    is 

K(R,pu)),G(R,pu))   =  ^.  (2.7) 

_  »  »     1^  *  _J 

From  the  definition  of  K,  Eq.  (2.5);  we  find  that 

K(R,pcj)  =  CO  -  p^/2m  -!-  Ze^/|R|  -  M(R,pu))  (2.8) 

where  M  is  the  transform  of  M(1,  2)  with  a   =0^^. 

The  equation  defining  M;,  Eq.  (2.2),  indicates  that  M  can  he  found  if 
both  G,  and  Gp  are  known.   But  knowledge  of  the  exact  G,  and  Gp  usually  means 
that  the  entire  many  body  problem  has  been  solved  exactly.   This  is  rarely 
feasible,  and  so  it  is  at  this  point,  specification  of  the  form  of  the  opera- 
tor M,  that  various  approximations  are  introduced.   For  example,  the  Hartree- 
Fock  approximation  vhich  we  used  in  I  results  from  the  choice 

M(l,2)  =  -i5(l-2)  /  d(3)v(l-3)G^(5;3^)  +  lv(l-2)G^(l,2"^). 

In  this  paper,  we  shall  choose  an  M  which  embodies  the  approximation 
used  by  Gell-Mann  and  Brueckner  in  their  study  of  the  correlation  energy  of 
the  high  density  electron  gas.   The  physical  content  of  their  theory  may  be 
stated  as  a  prescription  for  the  construction  of  the  two  particle  Green  s 
function.   The  resiilts  of  the  original  G.B.  paper  can  be  obtained  from  a  two 

Q 

particle   Green's   fixnction  of  the  form 

G2(1,2;1'2')   =   G^(1;1')G^(2,2')-G^(1,2')G^(2,1')  .      -. 

+  i/hyd(5)d(i|)G^(l;5)Gj3,l')Z^(5,'^)G^(2,4)Gjli,2') 

-   i/h^  d(5)d(U)G^(l;5)G^(3,2')  v(5-^)G^(2;i|)G^(i+;l')  (2.9) 


where  G  is  the  single  particle  Green's  function  in  the  absence  of  inter- 
action  and  Z/(l,2)  is  the  effective  interaction,,  introduced  by  Hubbard^ 
vhich  satisfies  the  integral  equation 

T/{1,2)    =  v(l-2)  -  i/4iyd(5)d(4)v(l-3)G^(5,4)G^(4, 5)^(^,2).       (2.10) 

The  last  term  in  Eq.  (2.9)  gives  rise  to  the  second  order  exchange 

(2) 
energy  (e^    in  G.B.)  which,  for  the  infinite  electron  gas,  is  independent 

of  density.   In  a  non-uniform  electron  gas,  this  energy  might  possibly  be  a 
function  of  density  gradient  and  higher  derivations  of  the  density.     The 
approximation  procedure  we  shall  be  using,  however,  treats  gradient  depen- 
dence as  a  higher  order  effect  than  density  dependence  and  in  the  low  order 

approximation,  with  which  we  shall  be  concerned,  only  the  density  dependence 

2 
of  e,   would  appear.   We  suspect,  then,  that  in  some  way  this  last  term  will 

contribute  a  constant  term  to  energies  we  may  encounter.   Let  us  suppress 

this  term  for  the  time  being;  it  will  become  (in  Section  V)  apparent  how  to 

allow  for  its  effect. 

It  will  be  possible  for  us  to  work  with  the  Green's  function  G,  rather 

than  the  "unperturbed"  Green's  function  G  .   Let  us  make  the  replacement  in 

(2.9)  and  (2.10).   Having  made  this  choice,  we  may  express  the  left  side  of 

Eq.  (2.2)  as      - 

rd(2)v(l-2)G2(l,2;i;2'*') 

=  rd(2)  J5(l-2)  rd(3)v(l-5)G^(5,5'^)  -  v(l-2)G^(l,2"') 

+  i/^  rd(3)d('+)G^(l;2)v(l-5)^^(2,i^)G^(3,^)G^(4;5'*')>  G^(2;l'). 

(2.11) 


Considered  as  a  function  of  the  difference  between  its  two  time  coordinates, 
the  single  particle  Green's  function  is  continuous  except  at  zero.   It  fol- 
lows, therefore,  that  the  symbol  +  as  a  superscript  of  a  numerical  argument 
of  a  given  single  particle  Green's  function  has  significance  only  if  the  two 
time  variables  of  this  given  function  are  tied  together  so  that  they  differ 
only  by  the  infinitesimal  which  the  symbol  denotes.   Thus,  in  the  last  term 
of  Eq.  (2.11)  the  factor  G  (l,2'*")G  (^,5)  could  replace  the  factor 
G  (1,2)G  (i+,5"^)  without  changing  the  value  of  the  integral.   Having  made  this 
replacement,  we  observe  that  the  multiplier  of  G  ( 1, 2  )  within  the  brackets 
is  just  the  right  side  of  Eq.  (2.10);  i.e.,  it  is  the  effective  interaction 
itself.   Finally,  comparison  of  this  result  with  Eq.  (2.2),  defining  the  self 
energy  operator,  reveals  that  we  have  chosen  an  M  defined  by 

M(l,2)  =-i6(l-2)  /"d(5)v(l-5)G^(5,5^)  +  i2^(l,2)G^(l,2").  (2,12) 

This  is  the  self  energy  operator  we  shall  use.   It  embodies  the  major  feature 
of  the  G.B.  theory  in  a  form  very  similar  to  its  Hartree-Fock  counterpart. 
In  Hartree-Fock  theory,  the  exchange  pcter;tial  expresses  the  change  from  the 
average  energy  brought  about  by  a  correlation  arising  solely  from  the  exclu- 
sion principle.   In  the  G.B.  theory  the  repulsion  arising  from  the  coulomb 
force  induces  a  correlation  which  augments  the  one  due  to  the  exclusion  prin- 
ciple.  Both  effects  are  included  by  substituting  the  effective  interaction 
for  the  coulomb  interaction  in  the  exchange  term.   The  transform  of  this  self 
energy  operator  is  readily  evaluated  and,  when  inserted  into  Eq.  (2.8),  gives 

n(r)e 


2  2  r  ^5,' 

K(R,pcu)    =  cu   -   p   /2jn  +  Ze    /|r|    -    2   /      — dr  ■...■••.     .,:.-..■; 

<«  ■«  J       |R-r|         **  '    '  ~ 

-    i(2jtn)         /    dp'dio'2^(R,p-p' ,o.'-,o')G(R,p'w')e  '^'^,    ..^     ,   .,  (2.15a) 
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.  n(r)  =  -iG^(r,t;r,t^).  (S.ljb) 

The  factor  of  two  in  the  third  term  on  the  right  of  Eq.  (2.13a)  arises  from 
the  spin  summation. 

Ill   EXPANSION  IN  POWERS  OF  ^ 

We  assumed  in  I  that  the  transform  of  the  single  particle  Green's  func- 
tion could  be  expanded  in  powers  of  ?i  as 

G(R,pu)  =*^  21!  ^"^  G  (R,pco)  (3.1) 

and  that  a  similar  expansion  for  the  space-time  form  of  the  Green's  function 

•V,-.   11 
IS  also  possible. 

G^(l,2)  =  2^  ^-^  G  (1,2)  (3.2a) 

.1  =  0 


p  dpdw 
G.(l,2)  =   /  -21-^  G  (R,pco)exp 


i(p.r-urt)/fi 


(3.2b) 


I   =  .^1-  ^2  *  =  ^1-  *2     S  =  ILV  i-g) 


We  shall  have  no  further  use  for  the  symbol  G  to  denote  the  single  particle 
Green's  function  without  interaction  or  for  the  symbol  G^   to  denote  the  two 
particle  Green's  function.   To  avoid  confusion  with  the  G.  used  as  expansion 

(J 

coefficients,  we  shall  delete  the  subscript  1  from  the  single  particle  Green's 
function  with  interaction,  and  shall  denote  this  latter  by  G(1,  2)  or,  some- 
times, simply  by  G.  "    ■   ' 

We  found  in  I  that  the  operator  0  and  the  kernel  K  could  also  be  expand- 
ed in  powers  of  -ft  as 


-  11 


ft^  9. 


(5- 5a) 


K  = 


li"^  K. 


(3.313) 


J=o 


These  expansions  may  be  inserted  into  Eq.  (2.7)  and  the  povers  of  h  separated 
in  such  a  manner  that  the  zeroth  order  approximation  is  non-vanishing.   The 
resulting  set  of  equations  will  be 


9  Ik  ,G  I  =  h"^ 


(5.^£ 


j+k+m=n 


"L  ~  ~ 


K^,G^ 


m 


=  0 


n  >  0. 


(3.^b) 


The  terms  0.  are  given  explicitly  in  the  Appendix  of  I.   The  coefficients  K. 
J  'J 

may  be  found  by  assuming  the  possibility  of  expanding  if   as  a  series  in 
powers  of  -h, 


?/(R,poj)  = 


*^  2/.(R,pu)), 


(5.5) 


inserting  this  series  into  Eq.  (2.15)^  and  separating  out  the  various  povers 
of  *.   We  find  tnat 


OOP    n   (r)e  dr 
K   (R,pu)    =  u)   -   p   /2in  +  Ze   /  R      -   2   / 

O    "V    w  "^  J 


R-r 


i(2rt)"      /   dp'dto'    y    (R,p-pSw-u)')G   (R,p'cj' )e^'^  ° 


o   .„  ■ 


o    «, 


(5.6a) 


p    n.(r)e   dr  , 

5.  =  -2  r  -^-^^^ — ^  -  -^ 


R-r 


(2rt)       -*      k+m=j 


dp'du)'  ?/'   G  e 


•     I    + 

~    ~     ico'o 

k  m 


(5.6b) 


n   (r)   =   -iG   (rt;r  t^) 

Jot  J      w       «i 

=   -i(27t)"      /    dpdcj  G.(R,pw)e'^*^° 


(5.6c; 
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■  '   A  calculation  of  the  coefficients  y  must  start  with  consideration  of 

n 

the  integral  equation  (2.10)  defining  the  effective  interaction.   Before 
starting  the  calculation,  however,  we  must  point  out  that  these  expansions 
in  powers  of  h  were  motivated  by  the  desire  to  include  the  effects  of  the 
uncertainty  principle  in  a  systematic  manner,  as  was  explained  in  I.   Al- 
though we  are  still  motivated  by  that  desire,  the  situation  here  is  somewhat 
more  complicated  than  that  in  I  because  of  the  appearance  of  h  in  a  way  that 
is  not  directly  related  to  the  uncertainty  principle.   We  refer  to  the  expli- 
cit appearance  of  h  in  equations  (2.9)  and  (2.10).   These  appearances  of  h 
have  their  origin  in  the  perturbation  expansion  and  diagrammatic  analysis 
upon  which  these  two  equations  are  based.   If  this  explicit  h  is  used  for 
expansion  purposes,  the  effective  interaction  will  be  decomposed  into  a 
series  of  terms  in  which  one  of  two  things  can  happen.   If  we  expand  so  that 
positive  powers  of  -ft  appear,  that  is,  so  that  the  series  is  dominated  by  its 
first  few  terms  as  ^   is  made  smaller,  then  the  lowest  order  term  which  ap- 
pears is  independent  of  the  strength  of  the  actual  force  between  particles 
and  there  is  no  way  of  recovering  the  exchange  term  in  that  limit  in  which 
correlation  effects  are  assumed  small.   On  the  other  hand,  if  we  expand  so 
that  negative  powers  of  h  can  appear,  then  the  result  is  a  series  containing 
arbitrarily  high  negative  powers  of ^,  which  looks  very  much  like  the  series 
of  perturbation  terms  which  G.B.  have  shown  how  to  sum.   Neither  of  these 
two  possibilities  is  at  all  in  accord  with  the  physical  interpretation  at- 
tached to  the  circumstance  that  the  effective  interaction  replaces  the  cou- 
lomb potential  in  the  exchange  term.   The  interpretation  of  this,  to  which 
we  called  attention  at  the  end  of  Section  II,  suggests  very  strongly  that 
the  effective  interaction  should  be  treated  as  an  entity  of  the  same  sort 
as  the  exchange  potential.   This  we  propose  to  do,  treating  the  effective 
interaction  here  in  the  same  way  as  we  treated  the  exchange  potential  in  I. 
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There,  however,  we  had  the  advantage  of  knowing  the  exchange  potential 
explicitly.   Here,  on  the  other  hand,  we  have  only  a  recipe,  the  integral 
equation,  for  calculating  the  effective  interaction.   This  recipe  calls  for 
that  h  inherited  from  perturbation  theory  as  one  of  the  ingredients.   This 
role  of  -^  should  be  suppressed  before  performing  any  expansion.   The  simplest 
method  of  suppressing  this  usage  is  to  replace  -fi  temporarily  by  the  symbol  X 
in  each  appearance  in  which  it  plays  the  role  of  the  perturbation  theory  ex- 
pansion parameter.   Instead  of  Eq.  (2.10),  for  example,  we  shall  have 


Y{1,2)    =  v(l-2)  -  /  d(^)Q(l,2)?/(4,2) 


(5.7) 


where 


;i,M  =  iA     d(5)v(i-5)G(5,MG(ii,5) 


(5.8) 


If  we  now  transform  these  equations  to  the  mixed  represent'ation,  the  trans- 
form of  the  integrals  may  be  expressed  in  terms  of  the  transforms  of  the 
factors  of  the  integrand  using  the  operator  0  once  again.   The  result  is 


z/'=  V  -  26 


Q  =  (i/x)9 


Q,v^ 


v,GG 


(5.9) 
(5.10) 


where 


GG  =   /  d(r^-r^)d(t^-tj^)G(5,^)G(4,5) 


X  exp  ■;  - 1 


I'^ly^k^-''^''^-^ 


/< 


(5.11) 


Again,  the  factor  of  two  in  Eq.  (5*9)  arises  from  the  spin  summation.   The 
integral  in  (5.11)  may  be  also  expressed  using  the  transform  of  the  Green's 
function  as 


Ih 


GG  =  (2jt^i)    /  dp'do)'  G(R,p'co')G(R,p+pS(d+u)') 


or,  using  the  expansion  of  the  transform,  as 


~    -2 

GG  =  ft 

1^      ^ 


_2  ^"^  "i^S'P^)  (3.12a) 


r>  dp'dw'    , 
"i  =   ^   /  -^^ r  ii^  G^(R,p'u)')  ft^  5  (R,IH-p',uH-to').       (5.12b) 

Lest  the  handling  of  the  powers  of  ii   here  seem  artificial,  let  us  note  that 

the  form  of  the  set  of  equations  (3.*+)  indicates  that  each  of  the  coefficients 

G,  contains  a  factor^  .   The  combination^  G.  which  appears  in  (5.12b)  has 
J  (J 

been  chosen  only  to  remove  this  h  dependence  from  the  terms  n.. 

J 

We  may  insert  the  expansion  of  ©  and  of  GG  into  Eq.  (5*10)  with  the  re- 
sult that 


Q(R,pw)  =    Zmj   '^    Q.(R,P«)  (5.15a) 

p6        ^ "  - 

Qi.   =  (kne^i/X)       2.   ®mC^/P^^'"n(S'r^^  '  (5.15b) 

Finally,  we  may  insert  the  expansions  of  V^©, and  Q  into  Eq.  (5.9)  and  sepa- 
rate the  varioxis  powers  of  -h.   The  coefficient  of  ^  then  satisfies  the  equa- 


tion 


^n  ^  2    2    e 

j+k+m=n   "^ 


=  (i^ne2/p2)5^^2-  (5.li^) 


Now  let  us  consider  each  of  the  equations  of  this  set  in  turn.   The 
n  =  0  equation  is 


1/    +2©  Q  ,1^ 
o     o|_^o  -^o 


=  0. 


■L} 


The  operator  0  acts  on  a  pair  of  transforms  to  produce  the  product  of 

these  transforms.   Hence,  the  above  equation  indicates  that  Z/  is  zero. 

o 

Every  term  in  Eq.  (j-l^)  vhich  contains  I/'    may  be  deleted  for  this  reason. 
Consider  now  the  n  =  1  equation  of  this  set.   Only  two  terms  survive  the 
deletion,  and  they  are  related  by  the  equation 


n  ^  29^ 


V^l 


=  0. 


Hence  7/     also  vanishes,  and  all  terms  containing  Z/^  may  be  deleted  from 
the  set  (5'1^)'   In  the  n  =  2  equation,  the  terms  surviving  may  be  re- 
arranged to  yield 

4ite  /p 

i?'(R,pw)  = .  (5.15) 

^  "  -     1  +  2Q  (R,pu)) 

The  higher  order  coefficients  i/.   may  be  expressed  in  terms  of  z/^  and  the 
other  Q  (n  <  j)  and  certain  differentiations  of  these  which  are  prescribed 
by  the  other  9  operators.  We  shall  not  need  these  higher  order  coeffi- 
cients here,  however.   This  completes  as  much  of  the  expansion  as  we  shall 
be  using  in  this  paper-   X  is  to  be  replaced  by  its  value  of  ft,  and  Q  may 

be  evaluated  using  (5.15b).   The  result,  Inserted  in  (5.15)",  gives 

2 

^2  =  ~ — ^^^^2 — r-  (5.16) 

p  +  Site  in  /a        ■■  ■■ 


IV  EQUATION  FOR  THE  CORRELATION  CORRECTION      ' 
Having  just  found  that  7f    and  j/"-,  vanish,  we  may  examine  Eq.  (5.6)  for 
K  and  K, .  '  '  ' 


..   K^  =  0,  -  p^/2m  +  Ze^/lRl  -  2e2  /  ^Ij:!-^         .  ,:  ^,,        (2,.ia) 


R-r 


16 


P   p  n  (r)dr 
K  =  -2e^   /   ^  "  "  .  (i+.lb) 

•^    |R-r| 


The  zeroth  and  first  order  equations  of  the  set  (5-^)  may  he  rewritten  as 

G  =^"^/k  (^.2a) 

o      '  o  ^ 

/k  .  (^.2b) 

o  ^ 


G,  =  -K,G  /k   -  qJk  ,G 
1     1  o'  o    1  o'  o 


These  equations  are  exactly  the  same  as  those  we  encountered  in  I  |j:.f.I;, 

Eqs.  (5-8),  (3.9),  and  (3.10)1.   We  found  there  that  2n  (r)  was  the  Thomas- 

— '  o  w 

Fermi  distribution  of  electrons^  and  that  K  and  G  both  vanished.   It  is 
convenient  to  express  (ij-.la)  and  ('+.2a)  in  the  following  form: 

G  (R,pa))  =*~^(cj  -  E)"""-  .  (^.3a; 

^3  V  ■» 

E  =  E(R,p)  =  p^/2m  +  ^^(R)  (^.3b; 


^^'  +   2e2 

p    n    ( r )  dr 
/         0    w      w 

R                    J 

R-r 

W\ 

\n\     m 

(R)  =  .££_+  2e^   /   °  "  "  .  {h.J>c] 


o 


d     is  then  the  usual  Thomas-Fermi  potential  for  the  atom. 
o 

The  integrations  which  yield  the  space  time  form  G  ,  and  the  density 
n  ,  as  indicated  by  Eq.  (5.2b)  and  (3-6c),  are  to  be  carried  out  as  if  w 
were  a  complex  variable,  with  a  path  of  integration  in  the  complex  w  plane 
which  is  displaced  from  the  real  w  axis  by  an  infinitesimal  amount.   The 
path  shown  in  Figure  1  lies  slightly  below  the  real  axis  for  -  00  <  u)  <  ^i, 
crosses  the  real  axis  at  cj  =  ja  and  lies  slightly  above  the  real  axis  for 
(a  <  CO  <  +  00.   This  choice  of  path  is  an  expression  of  the  fact  that  the 
system  is  in  the  ground  state,  as  explained  in  I.   The  constant  n  is  the 
chemical  potential,  approximately  the  energy  required  to  add  another  elec- 
tron to  the  system. 


I6c 
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Fig.  1.   Path  of  integration  in  the  complex  w  plane  for  recovering  G  from  G 


'  -■'!  ' . 


al 


(5.4; 


[  ^ 
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The  identity  of  these  equations  with  those  in  I  means  that  in  this 
case  too,  the  zeroth  order  result  is  the  Thomas -Fermi  distribution  and  that 
corrections  due  to  inhomogeneity  effects,  exchange,  and  correlations  are 
not  felt  until  the  second  order  in  -h.   Our  main  interest  is  to  obtain  an 
equation  for  the  first  non-vanishing  correction.   This  may  be  done  as  fol- 
lows:  Write  Eq.  (5-6b)  for  K  as 


K„(R,pu)  =  -</>J'R)    -   X(R,pu)) 


■T(R)  =  2e^ 


X(R,pw) 


n  (r)dr 


x|R-r| 


dp' 
(2«)- 


^  7/-(R,p-p',c.-u)')S  (R,p'co')e"'"'° 


(l+.l+a) 
(i+.l+b) 


/^    VVf   ryy 


We  may  rearrange  the  n  =  2  equation  of  the  set  (j.'+b)  so  as  to  exhibit  G^ 


Go  =  -K  G  /k  -  e„ 
2     o  o'  o    2 


K  ,G 
o  o 


/K. 


By  using  the  specific  fonn  for  K  and  G  given  by  Eqs.  (^.l)  and  (^.3),  this 
can  also  be  written 


^   Gg  =  Kg/Cw-E)^  +  (w-E)"^e2 


(w-E),(u)-E) 


1- 


(^.5) 


The  integral  of  Gp  indicated  by  (5-6c),  using  the  path  described  above, 
yields  the  density  np(R) . 

-,  n         n    A  •   +   K„(R,pw) 

(-2nh)^n  (R)  =  /  (iP  /   pTT  ^      7? 

2-    J   -.J   2rti         [u,-E(R,p)J^ 


dp 


du) 
2^ 


iu)0   ,   T-,N-1„ 
e     (u)-E)   e^ 


(w-E),(w-E) 


{^.6) 
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The  second  integral  in  Eq.  (^.6)  has  been  evaluated  in  I  [c.f.  I, 
Eqs.  (^.8),  (^.9),  (^+.10),  and  (4.11)].   It  is  equal  to 


12Pj,(R) 


d«^  \2 
o 

dR 


C+.T) 


where  p„(R)    is   the   Fermi   momentiim,    defined  by  the   equation 
F 


E[r,p    (R)]    =   n 


{k.d) 


Consider  now  the  first  integral  in  Eq.  (4.6).   The  substitution  of  Eq. 
(4.5)  and  (4.4)  in  this  expression  yields 


J   f-r.\     r  1   /   du)    iwo 
-^^W        dp  /  ^  e 


[L-E(R,p)] 


/■o  -K-N"?    Tj    J    ,     r    du)        dw'      iwo        iw'o      rw-E(R,p)'] 
+   (2.h)         /    dpdpj    /     -^     ^e  e  L       ^  J 


X   ru)'-E(R,p' )]""""  'i/jR,^'  ,^-^') 

^   »«   VM   W 


(^•9) 


The  use  of  the  special  path  for  the  u  and  w'  integrations  becomes  some- 
what inconvenient.   It  will  be  easier  for  us  to  integrate  u)  and  cj'  along 
their  real  axes,  provided  we  compensate  by  adding  small  imaginary  parts  to 
the  w-E  denominators  so  that  the  poles  lie  above  or  below  the  path  of  inte- 
gration in  the  correct  manner.   This  may  be  done  by  replacing  each  E(R,p)  in 
(4.9)  by 

A(p)=E(R,p)  -  ia£(p-pp)  ;   a  ^  O"^  (^-lO) 


e(x)  =  ±1  ;  X  J  0 


€(0)  =  0  . 


It  is  now  possible  to  introduce  new  variables 


T)  =  co-cj' 


(^.11) 
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into  the  second  Integral  in  (^.9)-  The  zeros  appearing  in  the  exponent  ^ 
are  to  be  interpreted  as  if  written  t  and  t' ,  each  of  which  approach  zero 
independently  through  positive  values.   The  second  integral  become 

(2.^)-5  Tdpdp'  f    |lL.|^e^[V(t.f)-.t']  ^^,^p.,,,,) 


;^.r  ;/'.,',;j.Q  eisuw 


-1 


X  W-A{v)Y      ^'-TACp')]'-"  .  -  ■•  -    L'-.'    (^.12) 

The  Ti  and  t]'    integrations  may  be  simplified  by  writing  the  factors  following 
V    in  the  form 


[t1-A(p)+A(p')]    \&'-ttA(p')]    -  &'-A(p)] 

-  lTi-A(p)+A(p')r^  h'-A(p)r^  .  (^-15) 

The  T]'  integration  associated  with  the  first  two  terms  of  (^.15)  may 
be  performed  by  contour  integration,  treating  r\   as  real.   In  the  limit 
t,t'  -  O"*"^  the  contribution  of  these  terms  to  the  expression  (■^.12)  will  be 

(2«fi)"5   Tdpdp'  r  1^  7„(R,P-p',ti)[t1-A(p)+A(p')]"^  -->..  e>n.-J- 


2rti  2' 


*»  JM  1^ 


7^°%(p-p')  -  e-^^°%(p  -p)" 


[y 


■J    a  ) 


where 


S(x)  =   1  X  Jo  . 


If  the  variables  p  and  p'  are  interchanged  in  the  second  term  of  the  inte- 
grand, and  if  the  substitution  of  -ri  for  t^  is  also  made  in  this  term,  the 
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result  "becomes 


(2«ft)-5      Tdpdp'    f^    |^e^^°     S(pp-p')rn-A(p)+A(p')]-^ 


1 
'    (R,p-p',Ti)    -         (R,p'-p,-Tl)     .  {^.IK 

*—        ^A       Wi        iVA  "^       y"JV\  1. 


We  shall  show  later  that  7/^  is  a  function  of  the  magnitude  |p-p'  |,  and 
is  an  even  function  of  r\.      Hence,  the  term  ahove  vanishes.   The  contribu- 
tion of  the  last  term  in  (^.13)  to  the  expression  (^.12)  may  be  combined 
with  the  first  term  in  (^.9)  after  the  integration  variable  in  this  latter 
expression  is  changed  from  co  to  ri'  .   We  obtain 


dp<;cz{   (R)+(2nh)-5    Tdp'     r^e-^^°^[T,.-A(p)+A(p')]"^i^(R.[p-p'].n) 


X 


ill   pin'°^r„..Ar^^i-2 

27t 


.     e^^  °    &'-A(p)]-^    .  (^.15) 


•Wii 


Again,  the  t|'  integration  may  be  performed  by  contour  integration,  com- 
pleting the  path  along  the  real  axis  by  the  infinite  semicircle  in  the  upper 
half  of  the  t)'  plane.  For  p  either  greater  than  or  less  than  p  ,  this  inte- 
gral vanishes.   For  p  =  p  ,  the  singularity  lies  astride  the  contour  and  the 

I* 

integral  becomes  infinite.   This  suggests  delta  function- like  behavior,  which 
leads  us  to  try  to  evaluate  the  integral  so  as  to  exhibit  this  delta  func- 
tion.  This  may  be  done  as  follows:  .  ,  • 


=  -  2  6(p_-p)  +  e(a)  (lt.l6) 

P         -f 
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The  term  of  order  a  arises  from  the  imaginary  addition  we  made  so  that  the 
real  axis  could  be  used  for  a  path  of  integration.   It  should  be  discarded. 
We  choose  spherical  polar  coordinates  for  p  and  p',  taking  the  polar  axis 
for  p  in  the  direction  of  p' .   The  integrand  of  (■^.15)  is  then  independent 
of  all  angles  except  0,  the  angle  between  p  and  p' .   Hence,  expression  ('+.15^ 
with  w  replacing  r\,    becomes 

Wmp  (R)^  (R)  +  [mp  (R)Ah^Jj(R)  (^.ITa) 


F'  ^^^2' 


"F 


J  =  /  p^dp  /  dcosG  /  ^  e"^'^"  [7o-A(pp)+A(p)]-^  "^^R,  [p^-p]  ,c.) .  (i|.17b; 
o 

We  return  to  Eq.  (^.6)  in  which  expressions  [k .Yj)  and  (^.7)  may  now 
replace  the  first  and  second  terms  on  the  right  hand  side.  It  is  possible 
to  eliminate  n  (r)  from  this  equation  using  Eq.  (4.ilb)  in  the  form 

The  resulting  equation  for  </)„{R)    is 

2  2 

V^  ^^(R)  -  ^  Pp(R)i/^2(R)  =  -f^     Pp(R)J(R) 
jth'  It  ^ 


2 
me 
+  — 


12rth5p^(R) 


.v^^  .^^"'° 


2~ 


o  2  \  dR  y 


(^.18) 


This  equation  determines  the  first  non-vanishing  correction  to  the  Thomas- 
Fermi  potential,  in  terms  of  the  Thomas-Fermi  potential  ^   (r),  the  Fermi 
momentiim  p„(R)  and  the  function  J(R),  which  arises  from  the  exchange  and 
correlation  effects.   In  contrast  to  the  situation  in  which  only  exchange 
is  considered,  it  will  no  longer  be  possible  to  express  the  dependence  of 
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J  on  p  in  any  simple  manner.   Still,  the  dependence  of  J  on  p   can  be 

r  r 

exhibited  much  more  explicitly  than  we  have  done  in  (^.17b)  and,  in  its 
more  explicit  form,  numerical  calculation  can  be  performed.   Accordingly, 
the  study  of  J(r)  is  the  subject  of  the  next  section. 

V  STUDY  OF  J(r) 

The  definition  of  J(R)  by  means  of  Eq.  (^.17b)  and  that  of  x^  by 
means  of  (5. lb)  leads  us  to  Eq.  (5'12b)  for  the  definition  of  H  : 


n  (R,pa))  =  (2n)"  *   /  dp'dw'G  (R  p' w' )S  (R.p+p' .t>H-co' )  (5-1) 


The  Green's  function  G  which  we  shall  need  in  order  to  evaluate  the 

o 

expression  above  is  only  partially  defined  by  Eq.  (4.5a).   The  missing 
part  of  the  definition  is  the  statement  that  the  path  of  integration  in 
the  complex  ui   plane  must  be  the  one  described  in  the  paragraph  immedi- 
ately following  Eq.  (4.5a).   This  fact  leads  to  manipulations  which  become 
quite  cumbersome  when  two  or  more  Green's  functions  are  present,  as  in  the 
expression  above.   We  may  avoid  these  manipulations  by  adding  small  imagi- 
nary parts  to  the  denominators,  i.e.,  by  again  replacing  the  E(R,p)  of  Eq. 
(4.5a)  by  the  A(p)  of  Eq.  (4.10)  and  using  the  real  co'  axis  for  the  path 
of  integration.   In  this  manner,  we  obtain 


(2rt)5in  =  -/"dp'  r|^  [7o'-E(R,p')+iae(p'-Pp)l"^ 


o    J   ~^  J     2ni 
lim 
a-*o 


ru)+w'-E(R,  Ip+p'  I  )  +  ioe (  1  p+p '  I -p  )"] 


f^      fM  1* 


The  w'  integration  may  be  performed  by  contour  integration,  with  the  result 
that 
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(2rt)\fi^  =    /      dp'|u-(p^+2p  p')/2m-2ia]"^ 
^1 


o 


J      dp"  [To-(p^+2p  p')/2m+2ia]""^       ^  (5-2) 


vhere   the   regions   R     and  R_   are   defined  by 
R-l  :         |p'  I    >  Pv     a'^^i      Ip+P'  I    <  Pt? 

Rp:   Ip'I  <  Pv  ^'^     Ip+p' I  >  Pt;,  • 

The  variable  of  integration  in  the  first  term  of  (5-2)  may  be  changed 
from  p'  to  q  =  -(p+p')  so  that 


(2rt)^iQ  = 
^      o 


-  [w  -  (p  +2p  q)/2m  +  2ia]'^>  .  (5-5) 


The  position  dependence  of  Q     comes  about  solely  because  the  Fermi  momen- 
tum is  position  dependent.   There  vill  therefore  be  a  marked  similarity  be- 
tween many  of  the  expressions  we  shall  use  and  analogous  expressions  which 

appear  in  the  study  of  the  uniform  electron  gas.   The  expression  (5- 3)  ap- 

12 
pears  in  the  work  of  DuBois,    for  example^  where  it  is  called  Q  ,  the  pro- 
pagator for  density  fluctuations  in  the  bare  pair  approxim.ation.   A  quantity 
proportional  to  it  is  denoted  '|y^(k;Lo)  in  Hubbard's  work.   The  mathematical 
formalism  we  shall  be  using  follows  closest  to  that  introduced  by  DuBois. 


-    2k 


If  we  use   dimensionless   variables   x,,    j,    and  u^    defined  by  x  =  p/p-rn^ 
y  =   f^/tPTT,  /2m)    and  u  =   q/p   ,    then 

(2rt)\fi^   =   2rtmp^(R)f(x,y)  (5Aa) 


i  du-  1  ^         •  1 


rt 


(5Ab) 


,  ... ,  2  2 

I     f-'  y  +  x  +  2x-u-    2ia        -y  +  x  +  2x-u  -    2ia 

I u+x I >  1 


The  form  of  f(x,y)  gives  us  some  valuable  information  even  before  the 
integration  is  performed.   We  note  that  f  behaves  as  l/y  as  y  -»■  oo;  that  it' 
is  an  even  function  of  y;  and  that  it  is  analytic  except  along  tvo  separate 
lines  in  the  complex  y  plane,  composed  of  those  points  for  which  one  or  the 
other  of  the  denominators  vanishes  for  some  u  within  the  domain  of  integra- 
tion.   These  lines  are  depicted  in  Figure  2.   The  function  f  depends  on  the 

\    .=-   ■  ■'■^    {."  --.  i.  Jt;-^rsrv.r!(^J^^rJ^;**v-V'>;v...!l;^^-*v^.-*.■^-■>r^!^-J■f.-Ii^.-'•■  .v  ■>.:•. ir.r-y  .,'^-;  jfirftK 

magnitude  of  x,  rather  than  on  its  direction. 

We  may  now  reconsider  the  expression  for  r/p  given  by  Eq.  (3>l6).   Using 
the  definition  (5.4),  we  have  i 

I   2   2/„ N         i 
h^e   p   (R)        I 

z/;(R,p  -p,co)  =  -^ ;  (5.5a) 

-^  ~      X  +  c(R)f(x,y)       \ 
where 

c(R)  =  2me^/(rthpp)  xi:JqfiO'.'   !.U-    (y..>c)1  ^o  3,-fuo  doxie'iS   .S;  -Sj:'¥    (5.5b) 

x  =  |p  -p|/p  (5.5c) 

y  =  co/(p//2m).  (5.5d) 

Those  assertions  about  ?/-  which  we  made  to  justify  setting  (4.1^)  equal  to 
zero  are  now  seen  to  follow  directly  from  some  of  those  properties  of  f(x,y) 
to  which  we  have  just  called  attention. 


-xU+2) 
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Fig.  2.   Branch  cuts  of  f(x,y)  in  complex  y  plane. 
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This  expression  for  xf    is  to  be  inserted  into  the  expression  (4.17b) 
defining  J(R) .   The  variables  of  integration  may  be  transformed  from  p, 
cos  9,  and  to  to  x,  y,  z,    where  z  =  |p|/p  so  that 


2        2 
X  =  1  +  z   -  2z  cos  9  . 


(5.6) 


Thus 


1+z 


J  = 


,   2 
4jte  p^ 


zdz 


xdx 


2rti 


-lyo 


1-z 


-1  r  2 


-1 


[y  -  1  +  z  -  iae(z-l)]"^  [x^  +  cf(x,y)]'^  . 


(5.7) 


The  order  of  integration  over  x  and  z  may  be  interchanged,  eind  the  z  inte- 
gration, whose  limits  become  1+x  and  |l-x|,,  may  be  performed.   The  result 
of  this  z  integration  will  depend  on  whether  or  not  the  point  z  =  1  lies 
within  the  range  of  integration.   Thus 


J  =  2ue\      r  xdx  r  g^  e-^y°  g(x,y)[x2  +  cf(x,y) 


-1 


(5.8a) 


.      y  +  ia        ,  y  -  ia  +  x(2+x)       ^  _ 
g  =  tn  — =^ TTT — r-  +  In  ^ :— i '-     :   x  <  2 


y  +  ia  -  x(2-x) 


y  -  la 


(5.3b) 


g  =  in 


y  -  ia  +  x(x+2) 
y  -  ia  +  x(x-2) 


X  >  2  , 


(5.8c) 


The  integral  over  y  can  be  considered  to  be  a  contour  integral  over  a 
path  enclosing  the  entire  lower  half  of  the  complex  y  plane.  Each  singu- 
larity of  the  integrand  will  yield  a  contribution  to  J(R)  and  so  these  singu- 
larities are  of  interest.   They  are  depicted  in  Figure  J-   I^gtire  5  indicates 

12  13 
the  so-called  plasmon  poles   '    contributed  by  the  second  factor  of  the  inte- 
grand, located  at  those  points  y  =  ±y  for  which 
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Pig.  5.   Singularities  of  the  integrand  of  Eq.  (5-8)  in  complex  y  plane. 
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X  +  cf(x,y^)  =  0.  (5.9) 

A  study  of  this  equation  indicates  that  for  values  of  x  and  c  for  which  the 
solution  lies  on  the  principal  branch  of  f^  that  solution  is  in  the  second 
and  fourth  quadrant  of  the  complex  y  plane.   We  go  into  this  matter  somewhat 
more  carefully  in  Appendix  A  of  this  paper. 

Instead  of  pursuing  the  steps  which  lead  to  the  conversion  of  the  y 
integral  to  integrations  over  the  finite  length  branch  cuts  in  the  lower 
half  plane,  we  shall  follow  a  course  which  will  show  the  connection  between 
this  integral  and  the  non  kinetic  part  of  the  separation  energy.   In  this 
way,  we  shall  be  able  to  avoid  the  necessity  of  a  numerical  evaluation  of 

that  integral  which  results  when  the  former  course  is  followed,  by  taking 

5  9  12  13 
advantage  of  n\imerical  results  others  have  obtained.  '    '      '  ye  rewrite 


(5.8a)  as 


2     r  dx   r  dy   ^-iyo 


'  -   ^-%    f    ST  ^    '      «(^'^' 


00 


+ 


n   2     /   dx   /  dy    -iyo    ,    v 


o 


cf(x.y) 


2 

X  +  cf(x,y) 


(5.10; 


The  first  integral  on  the  right  may  be  evaluated  by  contour  integration. 

The  only  singularity  of  the  integrand  in  the  lower  half  of  the  y  plane  is 

the  cut  which  extends  from  the  axis  to  x(2-x),  and  this,  only  for  x  <  2. 

(See.  Fig.  5)'   The  discontinuity  of  the  integrand  across  the  cut  is  -2iti, 

and  thus,  the  term  is  •        .     . 

2      x(2-x) 

J^^)  =  -2ne%p  f    f     f  ^^=   -''"^S  •  ^5.11) 

o      o   ■ 

In  the  second  term,  the  exponential  serves  to  inform  us  that  the 
singularities  in  the  lower  half  plane  contribute  to  the  y  integration  - 


i.e.,  that  the  integral  over  the  infinite  semicircle  on  the  lower  half  of 
the  y  plane  vanishes.   However,  since  f(x,y)  behaves  as  l/y  for  large  y 
(see  5.''+b),  this  exponential  is  not  needed,  and  can  be  set  equal  to  unity. 
The  path  of  integration  may  now  be  rotated  from  the  real  axis  to  the 
imaginary  axis  without  crossing  any  singularities.    Rotation  of  the  path 
of  integration  away  from  the  real  axis  gives  y  such  a  large  imaginary  part 
that  the  small  imaginary  additions  ia  which  appear  in  (5.8)  and  (5'^l3)  will 
have  no  effect.   When  a  is  set  equal  to  zero  in  Eq.  (5-^b)  defining  f(x,y), 
the  two  terms  of  the  integral  cancel  over  that  region  for  which  the  re- 
striction lu+xl  >  1  is  violated.   Hence,  this  restriction  may  be  dropped 
and  the  integration  performed,  yielding 


f(x,y)  =1+2^ 


-J 


x-f  ^ 


^n 


x(x+2)  +  y 


x(x-2)  +  y 


2x 


1?V^-  X 


In 


x(x+2)  -  y 
x-2)  -  y 


(5.12) 


This  form  is  closely  related  to  Lindliard's  dynamic  dielectric  constant 
€(k,Lo),   a  function  of  wave  number  k  and  frequency  cj.   In  fact, 

e(k,co)  =1  +  ^  f(x,y).  (5.15) 

We  may  now  insert  (5-15)  into  the  second  term  of  (5-10)  and,  on  setting 
the  exponential  there  equal  to  unity  and  rotating  the  path  of  integration 
from  the  real  to  the  imaginary  axis,  obtain 


r(2) 


2 


dx 

X 


— -  /   dy  g(xiy) 


;(l^u)) 


(5.1M 


We  note  that  Eq.  (5-12)  indicates  that  f(x,y*)  =  f*(x,y),  and  that  f(x,-y) 
=  f(x,y).   This  establishes  the  reality  of  f,  and  consequently  of  €, along 
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the  imaginary  axis.  We  also  note  that  setting  a  =  0  gives  (5. 8)  a  form 

«(^'^'=  '°^Tf[^  (5.15) 

•With  the  property  that  g(x,y*)  =  g*(x,y).   Hence  integration  of  the  product 
of  g(x,y)  with  an  even  function  of  y  along  the  imaginary  y  axis  will  not 
contain  the  contribution  of  the  imaginary  part  of  g(x,y),  and  we  may  replace 
g  by  its  real  part :  -      . 


Re  g(x,iy)=  -  In 


y  +  X  (x+2) 


2     ?     ? 
y^  +  x''(x-2)^ 


n  n  -  x/2 

/  du  .    (5.16) 

-^  n     f  Xv2    -  y  >2 


This  expression  may  be  inserted  into  (5.1^).   Introducing  v  =  y/2x,  z  =  x/2, 
we  have 

CM  00  1 

j(2)  ^   .l,e2p    r  ^^  f    d,  {j  -   1)    f    d^  %  "    .  (5.17) 

-00    o  _1     "-^^   ' 

We  have  achieved  the  desired  connection  between  our  quantity  J  and  a 

quantity  of  physical  significance  which  appears  in  the  study  of  the  infinite 

15 
electron  gas .  Quinn  and  Ferrell   compute  the  correlation  energy  of  an  in- 
finite electron  gas  by  studying  the  self  energy  of  an  electron  near  the  top 

of  the  Fermi  sea.   Comparison  with  their  work  reveals  that  J^   is  hn^ 

(2) 
times  the  exchange  energy  of  an  electron  at  the  Fermi  surface,  and  that  J^ 

is  h-n'^i   times  the  correlation  energy  of  an  electron  at  the  Fermi  surface. 

The  energy  required  to  remove  a  particle  of  momentum  p  is  called  the 

F 

separation  energy.   In  a  xiniform  system,  the  separation  energy  is  composed 
of  kinetic  energy,  potential  energy,  exchange  energy,  and  correlation  en- 
ergy -  that  is,  designating  the  separation  energy  by  -E  , 

D 

2 
Eo  =  P-o  /2m  +  E   (p„)  +  E    (p„). 
S    F  '      ex^-^F     corr^-^F 
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We  have   found  now  that 

J(R)    =  hT,\}    Eg   -   pp^/2m    )       .  (5.18) 

This  is  a  result  we  encountered  before  in  I,  although  we  did  not  recog- 
nize it  there.  The  equation  in  I,  (Eq.  ^+.15),  can  be  written  in  the  form  of 
Eq.  (^.l8)  of  this  paper  so  that  a  quantity  in  I  can  be  defined  which  corre- 
sponds to  the  J  of  this  paper.   This  quantity  turns  out  to  be  equal  to 

2  2 

-^ite  p  ,  which  is  ^jt  h  times  the  non  kinetic  part  of  the  separation  energy 
r 

for  an  infinite  electron  gas  in  the  Hartree-Fock  approximation. 

This  fact  suggests  the  way  to  include  the  effect  of  the  term  which-  we 

deleted  from  Eq.  (2.9)-   Let  us  recall  that  when  one  calculates  the  non 

kinetic  part  of  the  separation  energy  of  an  infinite  electron  gas  using  the 

theory  of  G.B.,  the  result  is  a  sum  of  terms,  each  term  arising  from  one  of 

the  members  of  the  right  hand  side  of  Eq.  (2.9).   It  seems  likely  that  the 

contribution  to  J  of  the  deleted  term  may  be  included  along  with  the  others 

merely  by  using  an  expression  for  E   in  which  the  contribution  of  this  term 

o 

appears.   Such  an  expression  can  be  derived  from  G.B.'s  foiroula  for  the 

1  o 
correlation  energy  per  particle.     It  also  seems  likely  that  the  way  to  go 

past  the  Gell-Mann  and  Brueckner  theory  is  to  use  an  expression  for  the 
separation  energy  which  has  validity  beyond  the  high  density  limit.   There- 
fore, we  shall  use  the  general  equation,  (5.18),  for  J(R)  without  specifying 

in  detail  how  the  separation  energy  is  to  be  calculated.   Inserting  this  in 

first 
Eq.  (4.18)  gives  the  equation  for  the^correlatlon  and  quantum  correction  to 

the  Thomas -Fermi  model: 


.=..r-4^,..^J.   ^^ 


2   1   ^5  i  ^F'^2      ,  5   )   S    2m 
^  Ti¥i  J  jth    ' 
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,,  „  2  ,     2m  /  '^o 


(5.19) 
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VI   SUMMARY  Affl)  DISCUSSION 

An  equation  has  been  derived  which  determines  the  first  correction  to 
the  Thomas-Fermi  potential  in  terms  of  the  Fermi  momentum  of  the  ordinary 
Thomas-Fermi  model  and  the  separation  energy  of  an  infinite  electron  gas. 
This  correction  describes  the  change  from  the  Thomas-Fermi  potential  due 
to  exchange,  inhomogeneity,  and  correlation  effects,  and  is  the  first  term 
of  a  series  vhose  higher  members  can,  in  principle,  be  found  using  the 

techniques  described  in  this  paper. 

7 
A  recent  paper  by  Lewis   describes  a  method  of  including  correlation 

and  exchange  effects  in  the  Thomas-Fermi  model.   His  approach  is  based  on 

the  assumption  that  the  electron  cloud  surrounding  that  atom  behaves  loca] 

ly  like  an  infinite  electron  gas,  and  therefore  that  a  local  separation 

energy  can  be  ascribed  to  each  point  in  the  atom.   This  separation  energy, 

he  reasons,  must  be  independent  of  position  if  the  electron  distribution 

is  to  be  equilibrium.   By  including  successively  the  contributions  of  the 

potential,  exchange,  and  correlation  energies  in  his  expression  for  the 

separation  energy,  he  obtains  successively  the  Thomas-Fermi,  the  Thomas- 

Ferml-Dirac,  and  the  Thomas-Ferml-Dirac-Correlation  model  of  the  atom. 

19 
Independently,  Tomishima   proposed  the  same  sequence  of  models  using 

a  slightly  different  approach.   This  time,  the  assumption  is  that  the  energy 
of  the  atom  is  given  by  the  volume  integral  of  an  energy  density  which  de- 
pends only  on  the  local  density  of  particles.   This  energy  is  then  minimized 
with  respect  to  variation  of  particle  density,  subject  to  particle  conserva- 
tion.  By  including  successively  the  various  contributions  to  the  energy 
density,  Tomishima  also  obtains  the  various  models  mentioned  before.   The 
relationship  between  his  work  and  that  of  Lewis,  to  which  it  is  equivalent, 
arises  because  of  the  way  the  separation  energy  is  related  to  the  energy 
density. 
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Both  authors  utilize  the  G.B.  expression  for  the  correlation  energy  at 
high  densities.   Both  recognize  that  the  expression  needs  modification  at 
low  densities  and  propose  changes  of  the  G.B.  expression  which  have  more  or 
less  reasonable  low  density  limits.   Some  such  prescription  would  have  to 
be  used  in  the  application  of  our  work.   Tomishima's  paper  gives  numerical 
results  which  show  that  the  equation  he  and  Lewis  propose  gives  an  electron 
distribution  rather  close  to  that  predicted  by  the  Thomas -Fermi-Dirac  equa- 
tion.  In  particular,  the  difficulty   of  the  discontinuity  in  density  at 
the  edge  of  the  atom  persists,  even  though  correlation  effects  are  included. 

A  virtue  of  the  teclinique  we  have  proposed  is  that  it  does  not  yields  in  the 

21 
no- correlation  case,  the  discontinuity  at  the  edge  of  the  atom.     It  is 

likely  that  the  distribution  will  behave  sensibly  at  the  edge  of  the  atom 

for  the  present  theoi'y,  in  which  the  effect  of  correlations  is  included. 

One  may  ask  how  our  technique  differs  from  that  of  Lewis  and  Tomishima 

for  those  more  or  less  extended  systems  in  which  inhomogeneity  would  not  be 

expected  to  play  a  major  role.   All  inhomogeneity  effects  can  be  suppressed 

in  our  theory  merely  by  discarding  the  0.  operators  for  j  >  0.   If  we  do 

J 

this,  we  can  sum  the  expansion  and  obtain  equations  almost  identical  to  those 
which  describe  an  infinite  medium,  the  sole  difference  being  that  the  Fermi 
momentum  now  depends  on  position.   Such  a  procedure  applied  to  the  technique 
in  I  does,  in  fact,  lead  to  the  Thomas -Fermi-Dirac  equation.   So  far,  how- 
ever, we  have  not  been  able  to  recover  the  starting  point  of  the  Lewis  and 
Tomishima  technique  by  ignoring  the  inhoraogeneities  and  suimning  the  expan- 
sion.  The  difficulty  seems  to  be  connected  with  the  fact  that  when  correla- 
tion effects  are  introduced,  the  electrons  no  longer  occupy  a  Fermi  distri- 
bution in  momentum  space.   It  may  be  possible  in  some  other  manner  to 
derive  the  starting  point  of  the  Lewis  and  Tomishima  theory  from  ours  by 
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ignoring  inhomogeneities .   We  suspect^,  however,  that  there  may  be  a  real 
difference  in  content  which  makes  it  impossible  to  do  so- 

We  should  like  to  acknowledge  the  stimulating  conversations  with 
Professor  Sidney  Borowitz  during  the  course  of  this  work. 
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APPENDIX  A 
THE  PLASMON  POLES 

The  integral  (5-^b)  defining  f(x^y)  can  be  evaluated  for  real  y  in  the 
limit  a  -*  0  by  splitting  the  integrand  into  its  real  and  imaginary  parts  us- 
ing the  relationship 


Lim  — ^—. — 
z  +  la 
a-*o 


The  result  is 


P  -  ±   iit5(z). 


Lim_^   f(x,y)  =  A(x,y)  +  i  Z(x,y)  (A. la) 

CCr+O 

y  real 

A  =  i  P     f         du  i  ^i +  ^ >  (A. lb) 

"^  ^  '^   y  +  x^  +  2x-u   -y  +  X  +  2x-u 

iu+xl  >  1 


du  ^    8(y  +  x^  +  2x-u)  +  6(-y  +  x^  +  2x-u);>  .    (A.Ic) 

lul  <  1     ^ 
I u+x 1  >  1 

12 

These  expressions  appear,  in  basically  the  same  form,  in  DuBois'  paper. 
He  credits  Ferrell  with  the  observation  that  tha  restriction  |u+x|  >  1  is 
superfluous  in  the  evaluation  of  A  because  the  two  terms  of  the  integrand 
cancel  over  the  region  ju+x]  <  1.  The  evaluation  of  each  term  of  E  in- 
volves  calculating  the  area  of  a  plane,  perpendicular  to  x,  whose  distance 
from  the  origin  depends  on  x  and  y  and  which  cuts  one,  both,  or  neither  of 
the  unit  sphere  which  bound  the  region  of  integration.   These  spheres  are 
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centered  on  the  origin  and  on  the  point  -x.  For  x  >  2,  the  spheres  fail 
to  intersect,  and  so  the  form  of  the  result  differs  from  that  for  Ixl  <  2. 
Upon  evaluating  the  integrals,  we  obtain  the  result 


A  =  i 

X 


X  + 


^  +  x^2 


tn 


y+xTx+g) 


y  +  x(x-2) 


-m-y 


In 


y  -  x(x+2) 

y  -  x(x-2) 


(A. 2a] 


Z(x  >  2;y)  =  0 


2x 


1  - 


Kf-^'^ 


;   x(x+2)  <  y 

;   x(x-2)  <  y  <  x(x+2) 

;   0  <  y  <  x(x-2) 


(A. 2b) 


e(x  <  2,y) 


■Si 

1 

1/ 

y 

2x 

"H 

X 

ny 

2x 

;   x(2+x)  <  y 


;   x(2-x)  <  y  <  x(2+x) 


;   0  <  y  <  x(2-x) 


e(x, -y) 


z(x,y). 


This  result  was  obtained  previously  by  Hubbard,   although  his  A  and  E 

2 
differ  from  ours  by  a  factor  involving  p„  and  x  .   The  utility  of  this  result 

r 

is  twofold.   Firstly,  integrals  over  y  involving  f  can  be  frequently  trans- 
formed to  contour  integrals  to  be  evaluated  by  shrinking  the  contour  down 
around  the  singularities  of  f .   Once  these  have  been  identified  the  limit 
a  -»■  0  may  be  taken.   It  will  turn  out  that  the  integrals  remaining  will  in- 
volve only  A  and  E   ,  to  be  integrated  over  finite  lengths  of  the  real  y  axis. 
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Secondly,  this  result  determines  the  appropriate  branch  of  the  function 
which  describes  f  when  a  5^  0,  and   y  is  anywhere  in  the  complex  plane.   That 
is,  the  integral  (5-^b)  for  f  can  be  evaluated  just  as  it  stands,  with  the 
result 


f(x  <  2,y)  =  -  /x-ia  +  ^ 


1 .  hi^ ,  ,f 


In 


y-2ia+x(2+x) 
y-2ia+x(2-x) 


i(y-2ia)tn 


y-2ia+x(2-x)" 
y  -  2ia 


+  i(y+2ia)ln 


y+2ia-x(2-x) 
y  +  2ia 


1  -  J(£f2  - .)' 


In 


y+2ia-x(2+x) 


y+2ia-x(2-x) 


(A. 5a) 


f(x  >2,y)  =^L 


2ia   1 
X    2 


1  -  |(^  .  y.f 


In 


y-2ia+x(x+2) 


y-2ia+x(x-2) 


1 
+  2 


"1  _  l(y±2ia  _  ^)2 
4^   X 


In 


y+2ia-x(x+2) 


y+2ia+x(x-2) 


(A. 5b) 


The  phase  of  the  logarithm  of  each  divisor  or  dividend  should  be  re- 
stricted to  -n  <  iZ^  <  +  rt,  and  (A. 5)  will  then  go  over  into  (A.2)  in  the 
limits  y  real  and  a  -*  0. 

This  form  (A. 3)  is  useful  in  locating  the  plasmon  pole  at  that  y^  for 
which 


X  +  cf(x,y^)  =  0 


Re  f(x,y^) 
Im  f(x,y^) 


2/ 
-X  /c 


=  0  . 


(A.i+) 
(A.5) 


V/e  regard  the  position  of  the  pole  as  being  determined  by  the  intersection 


of  the  curves  in  the  y  plane  defined  by  (A.4)  and  (A.5).   Since  the  curve 
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(a. 5)  is  independent  of  c,  it  is  convenient  to  study  this  curve,  and  to 
regard  the  position  of  the  pole  along  it  as  being  determined  by  (A.^).   The 
general  behavior  of  the  curve  is  depicted  in  Fig.  h-.      The  relevant  features 
of  the  curve  are  these: 

For  y  far  from  the  imaginary  axis,  the  curve  is  asymptotic  to  a  straight 
line  whose  angle  with  the  real  axis  is  proportional  to  a. 

For  y  closer  to  the  imaginaiy  axis,  the  curve  approaches  the  end  of  the 
branch  cut  of  f(x,y),  and,  crossing  the  cut,  continues  on  a  different  sheet 
of  the  function.   The  distance  of  this  crossing  from  the  end  of  the  cut  is 
also  proportional  to  a. 

The  plasmon  pole  on  the  sheet  of  the  function  which  is  of  interest  to 
us   always  lies  in  the  second  and  fourth  quadrant,  if  it  is  present  at  all, 
and  thus  the  rotation  of  the  axis  used  in  Section  V  can  be  performed  without 
crossing  any  singularities . 

In  the  limit  a  -^-  0,  the  plasmon  pole  in  the  fourth  quadrant  lies  just 

o 
below  the  real  axis  to  the  right  of  the  cut.   As  x  /c  is  varied,  this  pole 

reaches  the  end  of  the  cut.   The  value  of  x  for  which  this  occurs  is  deter- 
mined by  (A.^)  by  setting  y  =x(x+2).   This  value  of  x,  which  depends  on  c, 

12  15 
is  called  x  ,  '        the  cut  off  momentum,  or  the  critical  momentum,  since  for 
c 

greater  values  of  x,  the  plasmon  pole  does  not  appear  but  seems  to  be  imbedded 

in  the  continuum.   As  x  is  reduced  from  values  greater  than  x  to  values  less 

^  c 

than  x  ,  the  plasmon  state  appears  at  the  end  of  the  continuum  represented  by 
the  cut,  and   detaches  itself  from  this  continuum.   This  behavior  was  some- 
what mysterious  when  the  plasmon  state  was  regarded  as  one  of  the  states 
arising  from  the  continuum.   Now  that  we  have  traced  its  position  onto  the 
non  physical  sheet  for  x  >  x  ,  this  mystery  is  dispelled. 


yoB.  - 


Fig.    h.      The   curve  Itn  f(x,y)   =   0. 
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APPENDIX  B 
TERMS  IN  THE  OPERATOR  9 


e^lK^G 


=  K  G 


e  [K,  G 


1  .   '^  /  SK  5g    Sk   hG 


xyz 


X       X 


.gL' 


K,G 


h 


/L~^ 


'  1  ^  a^K  a^G 


h^   a^ 


G 


2  I  ;,  2   .2 
xyz    \  ox   op 


Sp   hx 


a^    5^G     _^    a^      a^G 

xyz   .^^^P^Py     ^Px^Py  ^^ 


a^      a^G 
axap     ayap 


The  full  0  is  given  by 


K,G 


Lim     exp 
R'->R 

P'-*P 


ill  y~'  (  a      a       a      a 
2    ^^=i  \  ^~   3p^  "  3p~   M^ 

xyz     \      X         X  X         X /_ 


X  K(R,pu))G(R',p'u)) 
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